Data Summary

exposure_name <- c("met_homovanillic_acid", "met_histidine", "met_hydroxyphenyllactic_acid")

pfas_c <- c("pfas_pfda", "pfas_pfhpa", "pfas_pfna", "pfas_pfuna", "pfas_pfoa")
pfas_s <- c("pfas_pfhps", "pfas_pfhxs", "pfas_pfos")

outcome <- c("mic", "daystomic", "daystomac", "rapid", "daystorapid","daystohyp")

table1::table1(~factor(mic) + factor(rapid), full_data %>% select(exposure_name, pfas_c, pfas_s, outcome))
Overall
(N=40)
factor(mic)
0 30 (75.0%)
1 10 (25.0%)
factor(rapid)
0 29 (72.5%)
1 11 (27.5%) 
table1::table1(~., full_data %>% select(exposure_name, pfas_c, pfas_s, outcome[c(2,3,5)]))
Overall
(N=40)
met_homovanillic_acid
Mean (SD) 0.100 (0)
Median [Min, Max] 0.100 [0.100, 0.100]
met_histidine
Mean (SD) 66.8 (7.42)
Median [Min, Max] 65.9 [51.4, 82.0]
met_hydroxyphenyllactic_acid
Mean (SD) 0.559 (0.503)
Median [Min, Max] 0.374 [0.100, 1.68]
pfas_pfda
Mean (SD) 0.252 (0.107)
Median [Min, Max] 0.240 [0.141, 0.530]
pfas_pfhpa
Mean (SD) 0.299 (0.247)
Median [Min, Max] 0.225 [0.0707, 1.10]
pfas_pfna
Mean (SD) 1.17 (0.645)
Median [Min, Max] 0.965 [0.290, 2.80]
pfas_pfuna
Mean (SD) 0.155 (0.128)
Median [Min, Max] 0.120 [0.0707, 0.790]
pfas_pfoa
Mean (SD) 3.16 (1.38)
Median [Min, Max] 2.90 [0.960, 9.30]
pfas_pfhps
Mean (SD) 0.273 (0.124)
Median [Min, Max] 0.245 [0.141, 0.710]
pfas_pfhxs
Mean (SD) 3.37 (3.04)
Median [Min, Max] 2.45 [0.890, 15.0]
pfas_pfos
Mean (SD) 9.13 (3.79)
Median [Min, Max] 8.20 [3.10, 22.0]
daystomic
Mean (SD) 3110 (1670)
Median [Min, Max] 3650 [298, 5510]
daystomac
Mean (SD) 3470 (1500)
Median [Min, Max] 4000 [629, 5510]
daystorapid
Mean (SD) 3130 (1560)
Median [Min, Max] 3290 [348, 5510]

Lucid Analysis

  • Exposure: met_homovanillic_acid,met_histidine,met_hydroxyphenyllactic_acid
  • Latent clusters:
    • Carboxylic acids: pfas_pfda, pfas_pfhpa, pfas_pfna, pfas_pfuna, pfas_pfoa
    • Sulfonic Acids: pfas_pfhps, pfas_pfhxs, pfas_pfos
  • Outcome: mic, daystomic, daystomac, rapid, daystorapid,daystohyp
  • Note: met_homovanillic_acid was not included in the analysis since unique value
G = full_data %>% select(exposure_name[-1]) %>% as.matrix() %>% scale()

Z =  list(Carboxylic = full_data %>%
            select(sample_id, all_of(pfas_c)) %>%
            column_to_rownames("sample_id") %>%
            as.matrix() |> scale() ,
          Sulfonic = full_data %>%
            select(sample_id, all_of(pfas_s)) %>%
            column_to_rownames("sample_id") %>%
            as.matrix() |> scale() )

Y = full_data %>% 
  dplyr::select(all_of(outcome)) %>%
  mutate_at(.vars = c(outcome),
            .funs = ~unlist(.) %>% as.numeric()) %>%
  mutate_at(.vars = c(outcome[c(2,3,5,6)]),
            .funs = ~scale(.))
  

df <- data.frame()
plot <- vector("list", length(length(outcome)))

for (i in c(2,3,5,6)){
  fit <- LUCIDusM::est_lucid(
    G = G, 
    Z = Z,
    Y = Y[,i],
    K = rep(2, length(Z)),
    family = "gaussian",
    max_itr = 200,
    useY = TRUE,
    modelNames = rep("EEV", length(Z)))
    fit_reordered <- reorder_lucidM(fit, reference = c(1,2))
    # fit_reordered <- fit
    est <- data.frame(
                     outcome = colnames(Y)[i],
                     Exposure1_to_Carboxylic =
                       fit_reordered$res_Beta$Beta[[1]][,2],
                     Exposure1_to_Sulfonic = 
                       fit_reordered$res_Beta$Beta[[2]][,2], 
                     Exposure2_to_Carboxylic =
                       fit_reordered$res_Beta$Beta[[1]][,3],
                     Exposure2_to_Sulfonic = 
                       fit_reordered$res_Beta$Beta[[2]][,3],
                     # Exposure3_to_Carboxylic =
                     #   fit_reordered$res_Beta$Beta[[1]][,4],
                     # Exposure3_to_Sulfonic = 
                     #   fit_reordered$res_Beta$Beta[[2]][,4],
                     Carboxylic_to_outcome =
                       fit_reordered$res_Delta$fit$coefficient[2],
                     Sulfonic_to_outcome = 
                       fit_reordered$res_Delta$fit$coefficient[3])
    df <- df %>% bind_rows(est)
    (plot[[i]] <- plot_lucid_in_parallel_plotly(fit_reordered,
                                      sankey_colors =
                                        sankey_colors,
                                      text_size = 12))
    names(plot)[i] <- paste0(colnames(Y)[i])
    htmlwidgets::saveWidget(plot[[i]], file = fs::path(dir_figure,str_c("lucid_in_parallel_", colnames(Y)[i],".html")))
}

for (i in c(1)){
  fit <- LUCIDusM::est_lucid(
    G = G, 
    Z = Z,
    Y = Y[,i],
    K = rep(2, length(Z)),
    family = "binomial",
    max_itr = 200,
    useY = TRUE,
    modelNames = rep("EEV", length(Z)))
    fit_reordered <- reorder_lucidM(fit, reference = c(1,1))
    # fit_reordered <- fit
    est <- data.frame(
                     outcome = colnames(Y)[i],
                     Exposure1_to_Carboxylic =
                       fit_reordered$res_Beta$Beta[[1]][,2],
                     Exposure1_to_Sulfonic = 
                       fit_reordered$res_Beta$Beta[[2]][,2], 
                     Exposure2_to_Carboxylic =
                       fit_reordered$res_Beta$Beta[[1]][,3],
                     Exposure2_to_Sulfonic = 
                       fit_reordered$res_Beta$Beta[[2]][,3],
                     # Exposure3_to_Carboxylic =
                     #   fit_reordered$res_Beta$Beta[[1]][,4],
                     # Exposure3_to_Sulfonic = 
                     #   fit_reordered$res_Beta$Beta[[2]][,4],
                     Carboxylic_to_outcome =
                       fit_reordered$res_Delta$fit$coefficient[2],
                     Sulfonic_to_outcome = 
                       fit_reordered$res_Delta$fit$coefficient[3])
    df <- df %>% bind_rows(est)
    (plot[[i]] <- plot_lucid_in_parallel_plotly(fit_reordered,
                                      sankey_colors =
                                        sankey_colors,
                                      text_size = 12))
    names(plot)[i] <- paste0(colnames(Y)[i])
    htmlwidgets::saveWidget(plot[[i]], file = fs::path(dir_figure,str_c("lucid_in_parallel_", colnames(Y)[i],".html")))
}

for (i in c(4)){
  fit <- LUCIDusM::est_lucid(
    G = G, 
    Z = Z,
    Y = Y[,i],
    K = rep(2, length(Z)),
    family = "binomial",
    max_itr = 200,
    useY = TRUE,
    modelNames = rep("EEV", length(Z)))
    fit_reordered <- reorder_lucidM(fit, reference = c(2,2))
    # fit_reordered <- fit
    est <- data.frame(
                     outcome = colnames(Y)[i],
                     Exposure1_to_Carboxylic =
                       fit_reordered$res_Beta$Beta[[1]][,2],
                     Exposure1_to_Sulfonic = 
                       fit_reordered$res_Beta$Beta[[2]][,2], 
                     Exposure2_to_Carboxylic =
                       fit_reordered$res_Beta$Beta[[1]][,3],
                     Exposure2_to_Sulfonic = 
                       fit_reordered$res_Beta$Beta[[2]][,3],
                     # Exposure3_to_Carboxylic =
                     #   fit_reordered$res_Beta$Beta[[1]][,4],
                     # Exposure3_to_Sulfonic = 
                     #   fit_reordered$res_Beta$Beta[[2]][,4],
                     Carboxylic_to_outcome =
                       fit_reordered$res_Delta$fit$coefficient[2],
                     Sulfonic_to_outcome = 
                       fit_reordered$res_Delta$fit$coefficient[3])
    df <- df %>% bind_rows(est)
    (plot[[i]] <- plot_lucid_in_parallel_plotly(fit_reordered,
                                      sankey_colors =
                                        sankey_colors,
                                      text_size = 12))
    names(plot)[i] <- paste0(colnames(Y)[i])
    htmlwidgets::saveWidget(plot[[i]], file = fs::path(dir_figure,str_c("lucid_in_parallel_", colnames(Y)[i],".html")))
}

writexl::write_xlsx(df, fs::path(dir_results, "LUCID/lucid_in_parellel_est.xlsx"))

Outcome: mic

plot[[1]]

Outcome: daystomic

plot[[2]]

Outcome: daystomac

plot[[3]]

Outcome: rapid

plot[[4]]

Outcome: daystorapid

plot[[5]]

Outcome: daystohyp

plot[[6]]

Outcome: Rapid. Exposure: 4 metabolites

G = full_data %>% select(exposure_name[-1], "met_glycine", "met_2_methyl_3_hydroxybutyric_acid") %>% as.matrix() %>% scale()

fit <- LUCIDusM::est_lucid(
    G = G, 
    Z = Z,
    Y = Y[,4],
    K = rep(2, length(Z)),
    family = "gaussian",
    max_itr = 200,
    useY = TRUE,
    modelNames = rep("EEV", length(Z)))

fit_reordered <- reorder_lucidM(fit, reference = c(2,2))

p <- plot_lucid_in_parallel_plotly(fit_reordered,
                                      sankey_colors =
                                        sankey_colors,
                                      text_size = 12)

htmlwidgets::saveWidget(p, file = fs::path(dir_figure,"lucid_in_parallel_rapid_with_4met.html"))
p